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Abstract 

Recently we have used a cellular automata model which describes the dynamics of a multi-connected network 
to reproduce the refractory behavior and aging effects obtained in immunization experiments performed with mice 
when subjected to multiple perturbations. In this paper we investigate the similarities between the aging dynamics 
observed in this multi-connected network and the one observed in glassy systems, by using the usual tools applied 
to analyze the latter. An interesting feature we show here, is that the model reproduces the biological aspects 
observed in the experiments during the long transient time it takes to reach the stationary state. Depending 
on the initial conditions, and without any perturbation, the system may reach one of a family of long-period 
attractors. The perturbations may drive the system from its natural attractor to other attractors of the same 
family. We discuss the different roles played by the small random perturbations ( "noise" ) and by the large periodic 
perturbations ( "immunizations" ) . 

1 Introduction 

In this paper we discuss a model for the evolution of the immune repertoire of B cells, which are responsible for 
the humoral immune responses. B cells belong to one of the main classes of white blood cells: the lymphocytes. 
These cells carry on their surface the order of 10 5 molecular receptors (proteins) and once activated they produce 
antibodies, which are copies of their molecular receptor. During the life time of an individual the immune system 
is able to produce the order of 10 11 different antibodies or different populations of B cells. The antigen (virus, 
bacteria, poison, cellular residue, etc) is not recognized as a whole but by its epitopes, which are patches on its 
structure that may be recognized by specific sites of the antibody molecules. By pattern recognition different 
antibodies will mark the epitopes of a given antigen, therefore forming a complex that will be eliminated by 
macrophages (another class of white blood cells) ^12- 

According to clonal selection theory |5] 0, elements that challenge the immune system will determine the 
populations (clones) of B cells that will proliferate: those populations will produce antibodies which will be able 
to recognize different epitopes of the specific antigen. The immune network theory however, is based on 

the fact that the antibodies (and molecular receptors) are able to recognize and to be recognized, and therefore 
during the immune response there are different types of interaction: antigen-antibodies and antibodies-B cells. 
In other words, when a given population of B cells is activated by the presence of a given antigen the produced 
antibodies will not only mark the specific antigens but also activate new B cell populations with complementary 
molecular receptors, which in turn will recognize them. The increase on the concentration of these complementary 
populations, on their turn, will maintain the proliferation of the antigen recognizing population, installing a 
feedback mechanism that will keep several populations activated. This kind of dynamics will generate a functional 
multi-connected network among different populations of B cells that will be dynamically regulated by mechanisms 
of activation and suppression. The network will then play an important role on the regulation of the immune 
responses. Although the immune network theory is part of the current immunological thinking, there are only few 
experiments supporting the interaction among clones with complementary receptors and the existence of such a 
network According to these experimental findings, if the network exists only 10 — 20% of the populations 

will belong to it, the rest of the immunocompetent populations remaining at rest. 

Recently we have successfully used a mathematical model 1161 [5] |7| |Sj (inspired in a previous one proposed by 
de Boer, Segel and Perelson |17|1 which takes into account the main features of Jerne's immune network theory, 
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to simulate experiments on immunization and aging performed with mice that could not be explained by the 
clonal selection theory . The simulations allowed to interpret the experimental results from the point of view of 
the immune network theory. 

The model allows to follow the evolution of the concentrations of the different populations of B cells in discrete 
shape space, a formalism which maps all possible molecular receptors of a given organism into points of a d- 
dimensional space. To each point (receptor) we associate a clone that corresponds to the population of B cells 
and antibodies characterized by this molecular receptor. The concentration of each clone will be described by 
a three-state automaton representing low, medium and high concentrations, and the interaction among different 
clones is based on complementarity. As far as we know, the model |17l 1161 1191 |H| |7j corresponds to the first 
successful attempt in describing the dynamics of the immune network as proposed by Jerne [!J| and recently it 
could be used to reproduce experimental results performed with mice §]■ Besides the biological implications of the 
results obtained in Ref. the dynamical behavior exhibited by this model in the biologically relevant parameter 
region is quite interesting by itself and should be better investigated. This region has been shown |Sj to exist 
for dimensions d > 2 and comprises a broad stripe near the transition between stable and chaotic behaviors, in 
which the model describes a multi-connected functional network In this region we found the majority of 

the populations in the resting state (low concentration) while the activated ones may reach 10 — 20% of the total 
number of populations. The activated populations are aggregated in clusters of different sizes, which fuse and 
split as time passes following an aggregation and disaggregation dynamics. Therefore the largest cluster at each 
time step is found in a different region \7\. 

In the experiments described in ref. 6 mice of the same linage are subjected to the immunization protocol as 
follows: the researchers inject ova by means of an intra-venal injection, wait for 14 days and measure the number 
of specific antibodies. Then they inject again the same amount of antigen, wait for 7 days and measure the amount 
of antibodies, and continue by repeating the same protocol every 7 days. In order to simulate the immunization 
protocol, the CA is subjected to specific perturbations by flipping chosen resting populations (low concentration) 
to their activated state (high concentration). Depending on the perturbation (damage) size, it may disappear 
after a few time steps or part of the damage (activated populations) may be incorporated to the network. In other 
words, the memory about the perturbations is due to the ability of the system to adapt (plasticity) and incorporate 
information about them. We have used two types of perturbations: random small ones, which correspond to the 
noise that mice are subjected to during the experiment, caused by the environment, and periodic large ones that 
will simulate the immunization protocol of multiple antigen perturbations under the same conditions After 
few presentations there is a saturation of the response of the system to the perturbation. During this process, the 
system incorporates new information and some of the previous information is lost, keeping the number of activated 
sites in the network almost unchanged. This kind of behavior was also observed experimentally in mice 10j , where 
saturation is related to a refractory behavior of the immune system. 

We have also observed aging effects on the dynamics of this system [S]. An older system is more rigid: the 
network loses plasticity to incorporate new information 1111 I12| . A recent study |13| has shown that the 
distribution of cluster sizes during the time evolution of the system has a characteristic cluster size (exponential 
behavior), but the distribution of persistence times (the period during which a given population remains activated 
and belongs to the network) exhibits a power law behavior. While the existence of a characteristic cluster size 
may be related to the loss of plasticity, the power law behavior of the persistence time may be associated to the 
memory generated by the dynamics of the system. 

The question we address here refers to how the learning process takes place dynamically and what is the 
cause of the loss of plasticity. The slow dynamics observed in this system presents analogies with the physical 
aging effects observed and reported on glass studies 1 141 115|: as the system gets older (ages) there is a loss of 
plasticity for structural or molecular relaxation and less changes are observed during the relaxation time. The 
mechanisms underlying the slow dynamics of glassy systems are the spatial (or geometric) disorder related with 
the difficulty to satisfy simultaneously all microscopic interactions, a characteristic called frustration. Glasses 
and spin glasses exhibit a rough energy landscape with many local minima which is responsible for slowing 
down the relaxation towards equilibrium. Their dynamics depends on the past history of the sample and under 
small perturbations they relax slowly. In our biologically motivated model the non-local interactions are based 
on complementarity (both perfect and slightly defective matches) and reflect the activation and suppression 
mechanisms in the complementary regions of the shape space. The analogue to frustration, in the network, is 
generated by the inability of the system to incorporate information and to satisfy the constraints of the activation 
and suppression mechanisms. We show that the CA dynamics gives rise to a large family of periodic attractors 
which are very robust, and could be regarded as the analogues of the minima in the glassy systems. The lost of 
plasticity and aging effects will therefore be related to the non-ergodicity in the phase space. 

In this paper the dynamics of the immunological responses in the network model are investigated using the 
common tools adopted in the study of glassy systems |T5 : . In particular we will focus on the relaxation of two- 
time autocorrelations of the B-cell populations, the structure of the attractors and effects of perturbations and 
immunizations. The paper is organized as follows: in section we describe the model and the procedures adopted 
in order to measure the correlations among different configurations of the system; in sectional we present and 
discuss the results obtained for small and large perturbations and in section [I] we present a summary, prospectives 
for future work and some conclusions. 



2 The Model 



The model under study here is a modified version |H||7||S| of the model proposed by Stauffer and Weisbuch |16|. 
which in turn was inspired in a previous model proposed by de Boer, Segel and Perelson |17| to describe the time 
evolution of the immune repertoire. It is a deterministic window cellular automaton model based on the shape 
space formalism 18] , which describes the interactions of B-lymphocytes and antibodies, and the main mechanism 
underlying these interactions, which is pattern recognition (lock- key interaction). The dynamics of the system 
is regulated by specific interactions involving complementary molecular receptors of the different clones. The 
memory about the relevant antigens, presented to the system during its past history, emerges from the dynamics 
of the system, rather than being stored in a static registry. 

To each point of a d-dimensional discrete lattice we associate a given receptor, which in turn will be described 
by d-coordinates representing important physical-chemical features of the receptor that differentiate one from the 
other |18|. Since clones are classified according to their molecular receptor, to each point r of the discrete shape 
space we associate a three-state automaton B(r, t) that will describe the concentration of its population over the 
time: low (B(f,t) = 0), intermediate (B(f,t) = 1) and high (B(f,t) = 2). 

The time evolution of the cellular automaton is based in a non-local rule: population B(f,t) at site r is 
influenced by the populations at site — r (its mirror image or complementary shape) and its nearest-neighbors 
(— r + 8r) (representing defective lock-key interactions). The influence on the population at site r due to its 
complementary populations is described by the field h(f,t): 

h(r,t) = Yl B ^'^ W 

where for a given r the sum runs over the complementary shape f" = — r and its nearest neighbors. Due to 
the finite number of states of the population B, the maximum value of the field h(r) is h max = 2(2d + 1). The 
updating rule is based on a window of activation, which describes the dose response function involved in B-cell 
activation ®E3[n|. There is a minimum field necessary to activate the proliferation of the receptor populations 
(di ), but for high doses of activation (greater than 62) the proliferation is suppressed. The updating rule may be 
summarized as: 

r,^,,)-! + ! if e 1 <h(r,t)<6 2 

( ' + > ~ \ B(r, t) - 1 otherwise 1 ' 

but no change is made if it would lead to B = — 1 or B = 3. We define the densities of sites in state i at time t 
as Bi(t) [i = 1, 2, 3). 

The initial configurations are randomly generated according to the following concentrations: -Bi(O) = £2(0) = 
a;/2, while the remaining L d (l — x) sites are initiated with B(f, 0) = 0. This model may exhibit stable or chaotic 
behaviors depending on the values of x, 9i and 62 — 9i- However it is on the transition region between the two 
behaviors that the model behaves like a multi-connected network [Hj. 

In order to simulate the immunization protocol performed in the mice experiments we have followed the 
procedure which is described in Ref. |5| and summarized below. We have adopted the scale of 5 time steps 
corresponding to 1 day While the system evolves according to the deterministic dynamics (eq. small and 
large perturbations can be produced, by setting the state of the chosen sites at B(r, i) = 2. 



Small Perturbations The small perturbations account for the immunological stimuli (noise) coming from 
the environment. The time interval between two consecutive small perturbations is a random number uniformly 
distributed between 1 and 100 time steps. Each perturbation corresponds to a random number of damages (from 
1 to 3) introduced at regions of resting populations (B — 0) which are randomly drawn (at every perturbation). 
The size of each damage may also vary randomly from 1 to 3 (onion-like) concentric layers around a central site 
(containing 7, 25 or 63 populations, respectively, in 3 dimensions). 



Large Perturbations The large perturbations correspond to the immunization protocol which starts at a 
predetermined age of the mice. Like in the experiments, we stimulate the system periodically (every 35 steps 
~ 1 week), and always in the same region (which is initially chosen at random but kept unchanged along the 
simulation). The damage size in this case corresponds to six layers (377 populations) around an specific site. 

Previous results on this model have shown [5] that the response to the immunizations presents a strong 
dependence on the initial time ("age") at which the periodic protocol starts (fitting experimental data extremely 
well for mice whose immunization protocol started at different ages). This has led us to the conjecture that the 
dynamical behavior of the system might be at least qualitatively similar to that of some glassy systems, despite 
the non-Hamiltonian nature of the CA dynamics. 

One of the quantities commonly used in the study of glassy systems is the two-time autocorrelation function 
between the system configurations at two given times t and t w . A common experiment in glassy systems consists 
in preparing the system at a high temperature and suddenly making a quench to a low temperature. Then the 
system is allowed to relax up to a waiting time t m , whose configuration is recorded. As the system continues 
to relax the autocorrelations between the instantaneous configurations at time t > t w and that at time t w are 
computed. The waiting time t w is called the age of the system in context of glassy systems. The monitoring of the 
two-time autocorrelations gives important insights on the relaxation process. The dynamics, whether stationary 



or not, can be readily recognized on the t w dependence of the autocorrelations, since in a stationary process 
two-time quantities depend only on time differences. Consequently by plotting correlations as a function of time 
difference t — t w for different waiting times it is possible to distinguish between an essential out of equilibrium 
process from a stationary one. The aging processes observed in glassy systems are then related to the lack of 
temporal invariance |15| . 

Inspired on this approach, we will analyze the multi-connected network dynamics by defining and analyzing 
quantities analogous to the two-time autocorrelations for the CA: 



C tot (t,t w ) = ±J25(B(?,t w ),B(?,t)) (3) 

r 

N 

C 22 (M») = B J w)N HB(r,t),2), (4) 

r|B(r,i ro ) = 2 

where N = L d is the total number of populations in the system. Ctot and C22 amount to normalized proximities 
(using Hamming distance as a measure) and from now on will be referred to simply as correlation functions. 
These quantities will be analyzed for different protocols: without any perturbation, with small perturbations 
(noise) and with large perturbations (immunization protocol) in order to differentiate the effects of the different 
kinds of perturbations. 

A complementary view of the long-term behavior of the system can be obtained by looking at the attractors 
to which the system evolves. In the present case it will consist mainly of cycles, which will be reflected by the 
periodicity of the correlations. For this purpose, we have obtained the return maps of the consecutive maxima of 
the correlation functions. Note that the period of the maxima thus obtained does not correspond to the period 
of the real system, which is at least twice as long |19|. 



3 Results 

3.1 Without perturbations 

In previous works |19ll51l7| it has been shown that depending on the initial concentration x of active populations the 
system may exhibit periodic or chaotic behavior. Systems with low initial concentrations of active sites (x < x c ) 
evolve to either fixed points or orbits with short periods, while for x > x c chaotic attractors appear. However, 
the biologically relevant region is in the transition region between these two behaviors, where the system reaches 
one of several periodic orbits (as will be seen below) with a very long period and after a long transient. From 
now on, all the results have been obtained using the same parameters adopted in Ref. d = 3, 0i = h max /3, 
82 = 2h m ax/3 and x — 0.26 (on the transition region). 

Without any perturbation, the system evolves after a transient time towards a cycle with a large period, 
as shown in the return map of Fig. Q We have also varied the waiting time from 10 to 100000 (not shown). 
The greater the waiting time, the closer to the attractor the system is, which is revealed by larger values of the 
autocorrelation functions. Once t m is larger than the transient, the time series for Ctot(t,tw) (and also C22{t,t w ) 
will include unity (see Fig.0 and will not change for larger values of t w . 

A typical result for the time evolution of the correlation functions for t w — 100 is shown in Fig. [5] The 
upper panel corresponds to the time evolution of the densities Bi (intermediate concentrations) and B2 (high 
concentrations), while the lower panel corresponds to Ctot{t, t w ) (open circles) and C22(i, t w ) (filled circles). Notice 
that while the concentrations relax to approximately constant values in a short time, Ctot and C22 take much 
longer to reach their attractors (note the logarithmic time scale) . In order to study how the system relaxes towards 
the attractor it is more convenient to make use of C22, since it measures the changes in the (more relevant) network 
of activated populations. In the case shown in Fig.|5|the transient time needed to attain the attractor is O(10 4 ) 
steps. It is necessary to point out, however, that this typical relaxation time is important only for the physical 
aspects of the dynamics. When mapped into the biological problem, it would correspond to ~ 5.5 years, which is 
much longer than the average life time of the mice used in the kind of experiment we simulated. Therefore the 
relevant behavior, from the biological point of view, happens to be in the transient of the model and not in its 
stationary state, a result which is interesting on its own. 

3.2 Random Small Perturbations 

How does the behavior of the system change when random small perturbations are produced on the parameter 
region used to simulate the real experiments performed with mice In order to investigate this issue, only the 
small perturbations described in Section [5] were produced, starting at time zero. 

In Fig.|3]we compare the time evolution of the densities and correlation functions obtained for single runs for the 
system without perturbations and with small random perturbations. Note that in the case with small perturbations 
the correlations initially follow those of the purely deterministic system, showing only small differences. After 
about 1 3 — 10 4 steps, however, they start decreasing faster, indicating some sort of cumulative effect that drives 
the system away from the region in phase space that it had approached until t w . These effects are more easily 
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Figure 2: Densities of intermediate and high concentrations vs. time (upper panel) and autocorrelations vs. t 
for L — 50 and t w = 100, without any perturbation. 



noted for C21. Moreover, the perturbations do not alter the behavior of the densities, as expected, since the 
number of activated populations is kept approximately constant by the self-regulatory mechanisms embedded in 
the dynamical rules. The changes are observed only in the populations that belong to the active network: in 
order to incorporate new information (new populations), older ones are deactivated. From the biological point of 
view, the realizations of the perturbations (small and large) differ from one individual to another, building the 
identity of each individual. At the end of life each individual will have a different history in terms of perturbations 
(antigen presentations) translated into the configuration of the populations belonging to the active network. This 
is nicely illustrated by Fig. 2] where two initially identical copies of a system undergo different realizations of the 
small perturbations according to the protocol described in section The Hamming distance between them grows 
on a long time scale, revealing the mechanisms behind the behavior of the correlations in Fig. |3] 

Returning to Fig. [3] it is important to stress that the changes observed in the concentrations for very long 
times (upper panel) are due to finite size effects. They are caused by the fact that all small perturbations are 
produced on regions of resting populations. For a finite system, after a long time all the possibilities will have 
already been explored. Increasing the size of the system, the changes on the densities disappear. This is shown 
in Fig. where we repeat the simulations, under the same conditions of Fig. |3J for a larger system (L = 100). 
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Figure 3: Densities vs. time (upper panel) and correlations vs. t — t w for L = 50 and t w = 100; without perturbation 
(circles) and with small perturbations (triangles). 
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Figure 4: Normalized Hamming distance between two initially identical configurations subjected to different sequences 
of small perturbations (three samples). 



In order to test out the robustness of these results, we have varied the parameters controlling the protocol of 
the random small perturbations. For instance, by increasing the maximum number of different perturbations from 
3 to 6 at each presentation, and/or by changing the maximum time interval between consecutive perturbations 
from 100 to 10, we observe the same qualitative behavior, with a faster decrease of the correlations — as expected, 
since in both cases the noise has been increased. 

What happens when the system's autocorrelations decrease due to the small random perturbations after it has 
reached a periodic attractor? Is the system driven to another cycle? To answer this question we have studied the 
stability of the cycles using the following procedure: we let the system evolve without perturbation towards its 
attractor for t w = 10000 time steps, after which we perturb the system with random small perturbations during a 
time interval Ati. Then we turn off the perturbations and allow the system to relax during another time interval 
At2 after which we obtain the return map of the correlations in the following 200 time steps. If we produce only 
one perturbation at t w (Ati = 1) the system remains in its original cycle, yielding a return map similar to that of 
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Figure 5: Densities vs. time (upper panel) and autocorrelations vs. t — t w for L = 100 and t w = 100 (the same 
conditions used in Fig. ; without perturbation (circles) and with small perturbations (triangles) . 
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Fig. Q Then, under the same conditions we repeat the simulations adopting now Ati = 1000 and At2 = 10000. 
We observe in Fig. |S| that the return map has changed slightly when compared to Fig. In particular, it no 
longer shows Ctot = 1 in the time series, but remains periodic, signaling that the system has shifted to a different 
cycle due to the perturbations. Apparently the periods observed in FigsQand|^]are the same. The distribution of 
periods and transient times are currently under investigation, results will be published elsewhere. Fig. |S| remains 
the same by increasing At2, which guarantees that the differences with respect to Fig.^are not a transient effect. 

According to these results, the role the noise plays, if allowed enough time to perturb the system significantly, 
is to drive the system from one attractor to a nearby one, which suggests that there is a family of periodic 
attractors which can be very close to one another in phase space. These effects, however, take place on a time 
scale which is longer [O(10 3 — 10 4 ) time steps] than the lifetime of the mice [O(10 2 — 10 3 ) time steps]. The small 
perturbations are therefore of little importance to the CA dynamics in the biologically relevant time scale, as 
suggested in previous work [5]. 

The behavior of the autocorrelation functions in Figs. |3]and|S]is reminiscent of what is observed in glassy 
systems. Here the CA approaches a periodic attractor, being thereafter deflected by the small perturbations. In 
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Figure 7: C22 (t + t w ; t w ) for different waiting times before small perturbations. From bottom to top, t w = 100, 5000, 



glassy systems, temperature drives the system from one local minimum to another, preventing it from getting 
stuck in local minima of the potential energy surface. As the physical age of the system (characterized by the 
value of the waiting time t w ) grows, it finds itself exploring deeper and deeper regions of a rough potential energy 
surface, diffusing towards equilibrium |20| . Due to the roughness of the potential energy, equilibrium is only 
attained on very long time scales and relaxation is very slow. The older the system is the more it gets confined 
to a restricted region of phase space and the time scales for relaxation get longer and longer. Eventually, if we 
wait long enough, the system equilibrates and the dynamics becomes stationary, losing sensitivity to the waiting 
time. This picture of an aging physical system is reminiscent to the loss of plasticity for adaptation in a living 
organism as it gets older. In either case, one observes a strong dependence on the waiting time t w , evident when 
measuring two-time quantities like autocorrelation functions and responses. Results for the CA model are shown 
in Fig. |7| where we see the decay of autocorrelations for three systems with different ages or waiting times. Note 
that the horizontal axis is the time difference between the total time and the waiting time. The three curves 
should collapse in the case that the dynamical evolution is stationary. 

3.3 Large immunizations 

What is the role of the large perturbations on the dynamics of the system? From previous studies we would expect 
the large perturbations to accelerate the aging process: while the random small perturbations would change the 
route to the natural attractor of the system, the large ones would reduce the transient time, a conjecture that 
would explain the loss of plasticity of the older mice UJ. The protocol adopted is the one described in Section [5] 
with six-layer perturbations every 35 time steps, always at the same sites. 

In Fig. |S] we compare the results obtained for the unperturbed system and those of the system subjected 
only to large immunizations starting at t w (note that the densities are now also plotted as functions of t — t w ). 
Somewhat surprisingly, the decrease of the correlations for large perturbations is small, when compared to the 
case of the small perturbations (compare with Fig.|3J. In hindsight, however, this can be understood because the 
large perturbations are always produced in the same region. For very long times the correlations will attain a 
stationary regime whose active populations contain at least part of the region involved in the immunizations UJ. In 
other words, the driving produced by the large periodic perturbations is of a completely different nature than that 
of the small perturbations. The large perturbations seem to play a selective role: the cycles that the system can 
reach are restricted to those that contain at least part of the populations incorporated during the immunization. 
According to this picture, the aging effects observed in Ref. 9 could be simply related to the exploration of the 
phase space: the older the system the less possibilities of choosing new cycles it would have. In Section 0] wc 
discuss this issue further. 

In Fig. |U] we compare the time evolution of the densities Bi and B2 and the autocorrelation functions for 
the system subjected only to large perturbations and for the system with both small perturbations and large 
immunizations. The results, as expected, confirm the dominance of the small perturbations, over the large ones, 
in driving the system faster to a different attractor. The increase of the densities around t ~ 10 4 for L = 50 
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Figure 8: Densities vs. time (upper panel) and autocorrelations vs. t — t w for a system with large immunizations, 
L = 50 and i w = 100. Triangles: without perturbations of any kind. Circles: both large immunizations and small 
perturbations. 
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Figure 9: Densities of active populations vs. time (upper panel) and autocorrelations vs. t — t w for t w — 100 and 
L = 50. Triangles: large immunizations only. Circles: both large immunizations and small perturbations. 



corresponds to the same finite size effects occurring in Fig. being associated to the active populations which 
have been incorporated into the network by means of the immunization protocol. 



4 Concluding Remarks 

We have analyzed the dynamics of a cellular automata model for the B-cell repertoire which has reproduced 
experimental results of immunizations on mice. Our analyses provide a broader context in which memory and 
plasticity take place, as discussed in Ref. 
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Figure 10: Hamming distance between 100 pairs of random configurations as a function of time (average and standard 
deviation). Inset: evolution for t > 100 shows that the HD attains a stationary regime for sufficiently long times 
(standard deviation not shown). 



Denning quantities analogous to the autocorrelation functions used in glassy systems, we have shown that 
the biologically relevant phenomena takes place in the transient regime of the model. The dependence of these 
functions on both t and t w (as opposed to the difference t— tw only) is reminiscent of the "aging" behavior observed 
in glassy systems, despite the fact that the underlying dynamics of both systems are controlled by completely 
different mechanisms. 

Starting from different initial conditions the deterministic dynamics takes the system towards a family of 
long-period attractors. When subjected to random small perturbations (Fig^J, the system is driven towards 
a new attractor of the family, revealing that most of the noise is assimilated. However, only part of the large 
perturbations is incorporated, due to the mechanisms of activation and suppression, leading to a saturation of the 
learning process 

From the biological point of view, the history of the mouse (sample) will be written by the different antigen 
presentations (random perturbations), starting from its "birth" (initial condition). Since the system is large but 
finite there is a maximum amount of information it may incorporate. The closer it is to its "destiny" attractor, 
the less information it is able to learn, since the deeper it already is in a given basin of attraction. Therefore the 
biological aging corroborated by the experimental results may be simply a consequence of this dynamical feature. 

From the results obtained up to now there are evidences that the purely deterministic dynamics is therefore 
non-ergodic. Despite their large number, however, the family of periodic attractors occupy only a fraction of 
the phase space. Evidences of this compression in phase space has been obtained in the following computer 
experiment: selecting randomly two initial configurations (with the same initial concentration), we measured the 
Hamming distance between them as a function of time as both systems evolve without any perturbation. In Fig. 1101 
we show the evolution of the average HD for 100 pairs (error bars are standard deviations). During the first 10 
time steps the HD decreases very rapidly, but for t > 100 we observe a quasi-stationary regime, reflecting the slow 
driving to the attractors (see inset). Note that this behavior is somewhat similar to that of the autocorrelation 
functions for small t w . 

Still focusing on the evidence of phase space compression, in Fig. llll we present the distribution of HD sampled 
from 500(500 — l)/2 pairs, for two consecutive time steps (t = 4000 and t = 4001). The distributions can be well 
described by a Gaussian, with a width that remains approximately constant (even for long times). Notice that the 
average value for t = 4001 is slightly larger than for t = 4000, reflecting the oscillatory behavior of the average HD 
as the pair of samples reaches their periodic attractors. It should be noted that, for large N, the Central Limit 
Theorem assures that randomly chosen initial configurations (with the same x) naturally give rise to a Gaussian 
distribution of the HD between any two of them, at time zero. Interestingly, the CA dynamics does not change 
the shape of the distribution, its only effect being to essentially shift the Gaussian towards lower mean values, on 
a long time scale. The spatio-temporal structure of the cycles, as well as their transients and basins of attraction, 
would be the object of further study and will be published elsewhere. 

This work was partially supported by CNPq, Faperj, Finep, Capes and Projeto Enxoval-UFPE. MC and RMZS 
acknowledge IF-UFF, where part of this work was done during their stay there. 
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Abstract 

Recently we have used a cellular automata model which describes the dynamics of a multi-connected network 
to reproduce the refractory behavior and aging effects obtained in immunization experiments performed with mice 
when subjected to multiple perturbations. In this paper we investigate the similarities between the aging dynamics 
observed in this multi-connected network and the one observed in glassy systems, by using the usual tools applied 
to analyze the latter. An interesting feature we show here, is that the model reproduces the biological aspects 
observed in the experiments during the long transient time it takes to reach the stationary state. Depending 
on the initial conditions, and without any perturbation, the system may reach one of a family of long-period 
attractors. The perturbations may drive the system from its natural attractor to other attractors of the same 
family. We discuss the different roles played by the small random perturbations ( "noise" ) and by the large periodic 
perturbations ( "immunizations" ) . 

1 Introduction 

In this paper we discuss a model for the evolution of the immune repertoire of B cells, which are responsible for 
the humoral immune responses. B cells belong to one of the main classes of white blood cells: the lymphocytes. 
These cells carry on their surface the order of 10 6 molecular receptors (proteins) and once activated they produce 
antibodies, which are copies of their molecular receptor. During the life time of an individual the immune system 
is able to produce the order of 10 11 different antibodies or different populations of B cells. The antigen (virus, 
bacteria, poison, cellular residue, etc) is not recognized as a whole but by its epitopes, which are patches on its 
structure that may be recognized by specific sites of the antibody molecules. By pattern recognition different 
antibodies will mark the epitopes of a given antigen, therefore forming a complex that will be eliminated by 
macrophages (another class of white blood cells) [1, 2]. 

According to clonal selection theory [2, 1], elements that challenge the immune system will determine the 
populations (clones) of B cells that will proliferate: those populations will produce antibodies which will be able 
to recognize different epitopes of the specific antigen. The immune network theory [3, 2], however, is based on 
the fact that the antibodies (and molecular receptors) are able to recognize and to be recognized, and therefore 
during the immune response there are different types of interaction: antigen-antibodies and antibodies-B cells. 
In other words, when a given population of B cells is activated by the presence of a given antigen the produced 
antibodies will not only mark the specific antigens but also activate new B cell populations with complementary 
molecular receptors, which in turn will recognize them. The increase on the concentration of these complementary 
populations, on their turn, will maintain the proliferation of the antigen recognizing population, installing a 
feedback mechanism that will keep several populations activated. This kind of dynamics will generate a functional 
multi-connected network among different populations of B cells that will be dynamically regulated by mechanisms 
of activation and suppression. The network will then play an important role on the regulation of the immune 
responses. Although the immune network theory is part of the current immunological thinking, there are only few 
experiments supporting the interaction among clones with complementary receptors and the existence of such a 
network [4, 5]. According to these experimental findings, if the network exists only 10 — 20% of the populations 
will belong to it, the rest of the immunocompetent populations remaining at rest. 

Recently we have successfully used a mathematical model [16, 6, 7, 8] (inspired in a previous one proposed by 
de Boer, Segel and Perelson [17]) which takes into account the main features of Jerne's immune network theory, 
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to simulate experiments on immunization and aging performed with mice [9] that could not be explained by the 
clonal selection theory . The simulations allowed to interpret the experimental results from the point of view of 
the immune network theory. 

The model allows to follow the evolution of the concentrations of the different populations of B cells in discrete 
shape space, a formalism which maps all possible molecular receptors of a given organism into points of a d- 
dimensional space. To each point (receptor) we associate a clone that corresponds to the population of B cells 
and antibodies characterized by this molecular receptor. The concentration of each clone will be described by 
a three-state automaton representing low, medium and high concentrations, and the interaction among different 
clones is based on complementarity. As far as we know, the model [17, 16, 19, 6, 7] corresponds to the first 
successful attempt in describing the dynamics of the immune network as proposed by Jerne [3] and recently it 
could be used to reproduce experimental results performed with mice [9]. Besides the biological implications of the 
results obtained in Ref. [9], the dynamical behavior exhibited by this model in the biologically relevant parameter 
region is quite interesting by itself and should be better investigated. This region has been shown [6] to exist 
for dimensions d > 2 and comprises a broad stripe near the transition between stable and chaotic behaviors, in 
which the model describes a multi-connected functional network [6, 7]. In this region we found the majority of 
the populations in the resting state (low concentration) while the activated ones may reach 10 — 20% of the total 
number of populations. The activated populations are aggregated in clusters of different sizes, which fuse and 
split as time passes following an aggregation and disaggregation dynamics. Therefore the largest cluster at each 
time step is found in a different region [7]. 

In the experiments described in ref. [9], 6 mice of the same linage are subjected to the immunization protocol as 
follows: the researchers inject ova by means of an intra-venal injection, wait for 14 days and measure the number 
of specific antibodies. Then they inject again the same amount of antigen, wait for 7 days and measure the amount 
of antibodies, and continue by repeating the same protocol every 7 days. In order to simulate the immunization 
protocol, the CA is subjected to specific perturbations by flipping chosen resting populations (low concentration) 
to their activated state (high concentration). Depending on the perturbation (damage) size, it may disappear 
after a few time steps or part of the damage (activated populations) may be incorporated to the network. In other 
words, the memory about the perturbations is due to the ability of the system to adapt (plasticity) and incorporate 
information about them. We have used two types of perturbations: random small ones, which correspond to the 
noise that mice are subjected to during the experiment, caused by the environment, and periodic large ones that 
will simulate the immunization protocol of multiple antigen perturbations under the same conditions [9]. After 
few presentations there is a saturation of the response of the system to the perturbation. During this process, the 
system incorporates new information and some of the previous information is lost, keeping the number of activated 
sites in the network almost unchanged. This kind of behavior was also observed experimentally in mice [10], where 
saturation is related to a refractory behavior of the immune system. 

We have also observed aging effects on the dynamics of this system [9]. An older system is more rigid: the 
network loses plasticity to incorporate new information [9, 11, 12]. A recent study [13] has shown that the 
distribution of cluster sizes during the time evolution of the system has a characteristic cluster size (exponential 
behavior), but the distribution of persistence times (the period during which a given population remains activated 
and belongs to the network) exhibits a power law behavior. While the existence of a characteristic cluster size 
may be related to the loss of plasticity, the power law behavior of the persistence time may be associated to the 
memory generated by the dynamics of the system. 

The question we address here refers to how the learning process takes place dynamically and what is the 
cause of the loss of plasticity. The slow dynamics observed in this system presents analogies with the physical 
aging effects observed and reported on glass studies [14, 15]: as the system gets older (ages) there is a loss of 
plasticity for structural or molecular relaxation and less changes are observed during the relaxation time. The 
mechanisms underlying the slow dynamics of glassy systems are the spatial (or geometric) disorder related with 
the difficulty to satisfy simultaneously all microscopic interactions, a characteristic called frustration. Glasses 
and spin glasses exhibit a rough energy landscape with many local minima which is responsible for slowing 
down the relaxation towards equilibrium. Their dynamics depends on the past history of the sample and under 
small perturbations they relax slowly. In our biologically motivated model the non-local interactions are based 
on complementarity (both perfect and slightly defective matches) and reflect the activation and suppression 
mechanisms in the complementary regions of the shape space. The analogue to frustration, in the network, is 
generated by the inability of the system to incorporate information and to satisfy the constraints of the activation 
and suppression mechanisms. We show that the CA dynamics gives rise to a large family of periodic attractors 
which are very robust, and could be regarded as the analogues of the minima in the glassy systems. The lost of 
plasticity and aging effects will therefore be related to the non-ergodicity in the phase space. 

In this paper the dynamics of the immunological responses in the network model are investigated using the 
common tools adopted in the study of glassy systems [15]. In particular we will focus on the relaxation of two- 
time autocorrelations of the B-cell populations, the structure of the attractors and effects of perturbations and 
immunizations. The paper is organized as follows: in section 2 we describe the model and the procedures adopted 
in order to measure the correlations among different configurations of the system; in section 3 we present and 
discuss the results obtained for small and large perturbations and in section 4 we present a summary, prospectives 
for future work and some conclusions. 
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2 The Model 



The model under study here is a modified version [6, 7, 8] of the model proposed by Stauffer and Weisbuch [16], 
which in turn was inspired in a previous model proposed by de Boer, Segel and Perelson [17] to describe the time 
evolution of the immune repertoire. It is a deterministic window cellular automaton model based on the shape 
space formalism [18], which describes the interactions of B-lymphocytes and antibodies, and the main mechanism 
underlying these interactions, which is pattern recognition (lock-key interaction). The dynamics of the system 
is regulated by specific interactions involving complementary molecular receptors of the different clones. The 
memory about the relevant antigens, presented to the system during its past history, emerges from the dynamics 
of the system, rather than being stored in a static registry. 

To each point of a d-dimensional discrete lattice we associate a given receptor, which in turn will be described 
by d-coordinates representing important physical-chemical features of the receptor that differentiate one from the 
other [18]. Since clones are classified according to their molecular receptor, to each point r of the discrete shape 
space we associate a three-state automaton B(r,t) that will describe the concentration of its population over the 
time: low {B{r,t) = 0), intermediate {B{r,t) = 1) and high {B{r,t) = 2). 

The time evolution of the cellular automaton is based in a non-local rule: population B(r,i) at site r is 
influenced by the populations at site —r (its mirror image or complementary shape) and its nearest-neighbors 
(— f + Sr) (representing defective lock- key interactions). The influence on the population at site r due to its 
complementary populations is described by the field h(f,t): 

h(r,t)= B ^',t) (1) 

r'£( — r+5r) 

where for a given r the sum runs over the complementary shape r' = —r and its nearest neighbors. Due to 
the finite number of states of the population B, the maximum value of the field h(f) is h max = 2(2d + 1). The 
updating rule is based on a window of activation, which describes the dose response function involved in B-cell 
activation [6, 16, 17]. There is a minimum field necessary to activate the proliferation of the receptor populations 
(0i), but for high doses of activation (greater than 82) the proliferation is suppressed. The updating rule may be 
summarized as: 

but no change is made if it would lead to B = — 1 or B = 3. We define the densities of sites in state i at time t 
as Bi(t) (i = 1, 2, 3). 

The initial configurations are randomly generated according to the following concentrations: Bi(0) = B 2 {0) = 
x/2, while the remaining L d (l — x) sites are initiated with B(r, 0) = 0. This model may exhibit stable or chaotic 
behaviors depending on the values of x, 81 and 82 —81. However it is on the transition region between the two 
behaviors that the model behaves like a multi-connected network [6]. 

In order to simulate the immunization protocol performed in the mice experiments we have followed the 
procedure which is described in Ref. [9] and summarized below. We have adopted the scale of 5 time steps 
corresponding to 1 day [9]. While the system evolves according to the deterministic dynamics (eq. 2), small and 
large perturbations can be produced, by setting the state of the chosen sites at B(f, i) = 2. 

Small Perturbations The small perturbations account for the immunological stimuli (noise) coming from 
the environment. The time interval between two consecutive small perturbations is a random number uniformly 
distributed between 1 and 100 time steps. Each perturbation corresponds to a random number of damages (from 
1 to 3) introduced at regions of resting populations (B = 0) which are randomly drawn (at every perturbation). 
The size of each damage may also vary randomly from 1 to 3 (onion-like) concentric layers around a central site 
(containing 7, 25 or 63 populations, respectively, in 3 dimensions). 



Large Perturbations The large perturbations correspond to the immunization protocol which starts at a 
predetermined age of the mice. Like in the experiments, we stimulate the system periodically (every 35 steps 
~ 1 week), and always in the same region (which is initially chosen at random but kept unchanged along the 
simulation). The damage size in this case corresponds to six layers (377 populations) around an specific site. 

Previous results on this model have shown [9] that the response to the immunizations presents a strong 
dependence on the initial time ( "age" ) at which the periodic protocol starts (fitting experimental data extremely 
well for mice whose immunization protocol started at different ages). This has led us to the conjecture that the 
dynamical behavior of the system might be at least qualitatively similar to that of some glassy systems, despite 
the non-Hamiltonian nature of the CA dynamics. 

One of the quantities commonly used in the study of glassy systems is the two-time autocorrelation function 
between the system configurations at two given times t and t w . A common experiment in glassy systems consists 
in preparing the system at a high temperature and suddenly making a quench to a low temperature. Then the 
system is allowed to relax up to a waiting time t w , whose configuration is recorded. As the system continues 
to relax the autocorrelations between the instantaneous configurations at time t > t m and that at time t w are 
computed. The waiting time t w is called the age of the system in context of glassy systems. The monitoring of the 
two-time autocorrelations gives important insights on the relaxation process. The dynamics, whether stationary 
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or not, can be readily recognized on the t w dependence of the autocorrelations, since in a stationary process 
two-time quantities depend only on time differences. Consequently by plotting correlations as a function of time 
difference t — t w for different waiting times it is possible to distinguish between an essential out of equilibrium 
process from a stationary one. The aging processes observed in glassy systems are then related to the lack of 
temporal invariance [15]. 

Inspired on this approach, we will analyze the multi-connected network dynamics by defining and analyzing 
quantities analogous to the two-time autocorrelations for the CA: 



where N = L d is the total number of populations in the system. Ctot and C22 amount to normalized proximities 
(using Hamming distance as a measure) and from now on will be referred to simply as correlation functions. 
These quantities will be analyzed for different protocols: without any perturbation, with small perturbations 
(noise) and with large perturbations (immunization protocol) in order to differentiate the effects of the different 
kinds of perturbations. 

A complementary view of the long-term behavior of the system can be obtained by looking at the attractors 
to which the system evolves. In the present case it will consist mainly of cycles, which will be reflected by the 
periodicity of the correlations. For this purpose, we have obtained the return maps of the consecutive maxima of 
the correlation functions. Note that the period of the maxima thus obtained does not correspond to the period 
of the real system, which is at least twice as long [19]. 

3 Results 

3.1 Without perturbations 

In previous works [19, 6, 7] it has been shown that depending on the initial concentration x of active populations the 
system may exhibit periodic or chaotic behavior. Systems with low initial concentrations of active sites (x < x c ) 
evolve to either fixed points or orbits with short periods, while for x > x c chaotic attractors appear. However, 
the biologically relevant region is in the transition region between these two behaviors, where the system reaches 
one of several periodic orbits (as will be seen below) with a very long period and after a long transient. From 
now on, all the results have been obtained using the same parameters adopted in Ref. [9]: d = 3, 81 = hm ax /3, 
82 = Ihmaxj'i and x = 0.26 (on the transition region). 

Without any perturbation, the system evolves after a transient time towards a cycle with a large period, 
as shown in the return map of Fig. 1. We have also varied the waiting time from 10 to 100000 (not shown). 
The greater the waiting time, the closer to the attractor the system is, which is revealed by larger values of the 
autocorrelation functions. Once t w is larger than the transient, the time series for Ctot{t,t w ) (and also C f 22(t,t w ) 
will include unity (see Fig. 1) and will not change for larger values of t w . 

A typical result for the time evolution of the correlation functions for t w = 100 is shown in Fig. 2. The 
upper panel corresponds to the time evolution of the densities Bi (intermediate concentrations) and B2 (high 
concentrations), while the lower panel corresponds to Ctotit, tw) (open circles) and (^(i, t w ) (filled circles). Notice 
that while the concentrations relax to approximately constant values in a short time, Ctot and C22 take much 
longer to reach their attractors (note the logarithmic time scale). In order to study how the system relaxes towards 
the attractor it is more convenient to make use of C22, since it measures the changes in the (more relevant) network 
of activated populations. In the case shown in Fig. 2 the transient time needed to attain the attractor is C(10 4 ) 
steps. It is necessary to point out, however, that this typical relaxation time is important only for the physical 
aspects of the dynamics. When mapped into the biological problem, it would correspond to ~ 5.5 years, which is 
much longer than the average life time of the mice used in the kind of experiment we simulated. Therefore the 
relevant behavior, from the biological point of view, happens to be in the transient of the model and not in its 
stationary state, a result which is interesting on its own. 

3.2 Random Small Perturbations 

How does the behavior of the system change when random small perturbations are produced on the parameter 
region used to simulate the real experiments performed with mice [9]? In order to investigate this issue, only the 
small perturbations described in Section 2 were produced, starting at time zero. 

In Fig. 3 we compare the time evolution of the densities and correlation functions obtained for single runs for the 
system without perturbations and with small random perturbations. Note that in the case with small perturbations 
the correlations initially follow those of the purely deterministic system, showing only small differences. After 
about 10 3 — 10 4 steps, however, they start decreasing faster, indicating some sort of cumulative effect that drives 
the system away from the region in phase space that it had approached until t m . These effects are more easily 
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Figure 1: Return map for maxima of the total correlation function C to t{t) for L = 50 and t w = 10000 (single run). 




noted for C 2 2- Moreover, the perturbations do not alter the behavior of the densities, as expected, since the 
number of activated populations is kept approximately constant by the self-regulatory mechanisms embedded in 
the dynamical rules. The changes are observed only in the populations that belong to the active network: in 
order to incorporate new information (new populations), older ones are deactivated. Prom the biological point of 
view, the realizations of the perturbations (small and large) differ from one individual to another, building the 
identity of each individual. At the end of life each individual will have a different history in terms of perturbations 
(antigen presentations) translated into the configuration of the populations belonging to the active network. This 
is nicely illustrated by Fig. 4, where two initially identical copies of a system undergo different realizations of the 
small perturbations according to the protocol described in section 2. The Hamming distance between them grows 
on a long time scale, revealing the mechanisms behind the behavior of the correlations in Fig. 3. 

Returning to Fig. 3, it is important to stress that the changes observed in the concentrations for very long 
times (upper panel) are due to finite size effects. They are caused by the fact that all small perturbations are 
produced on regions of resting populations. For a finite system, after a long time all the possibilities will have 
already been explored. Increasing the size of the system, the changes on the densities disappear. This is shown 
in Fig. 5, where we repeat the simulations, under the same conditions of Fig. 3, for a larger system (L = 100). 
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Figure 3: Densities vs. time (upper panel) and correlations vs. t — t w for L = 50 and t w = 100; without perturbation 
(circles) and with small perturbations (triangles). 
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Figure 4: Normalized Hamming distance between two initially identical configurations subjected to different sequences 
of small perturbations (three samples). 



In order to test out the robustness of these results, we have varied the parameters controlling the protocol of 
the random small perturbations. For instance, by increasing the maximum number of different perturbations from 
3 to 6 at each presentation, and/or by changing the maximum time interval between consecutive perturbations 
from 100 to 10, we observe the same qualitative behavior, with a faster decrease of the correlations — as expected, 
since in both cases the noise has been increased. 

What happens when the system's autocorrelations decrease due to the small random perturbations after it has 
reached a periodic attractor? Is the system driven to another cycle? To answer this question we have studied the 
stability of the cycles using the following procedure: we let the system evolve without perturbation towards its 
attractor for t m = 10000 time steps, after which we perturb the system with random small perturbations during a 
time interval Ati. Then we turn off the perturbations and allow the system to relax during another time interval 
At 2 after which we obtain the return map of the correlations in the following 200 time steps. If we produce only 
one perturbation at t w (Ati = 1) the system remains in its original cycle, yielding a return map similar to that of 
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Figure 5: Densities vs. time (upper panel) and autocorrelations vs. t — t w for L = 100 and t w = 100 (the same 
conditions used in Fig. 3); without perturbation (circles) and with small perturbations (triangles). 
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Figure 6: Return map for maxima of total correlation function, for L = 50, t w = 10000, Ati = 1000 and At2 = 10000. 



Fig. 1. Then, under the same conditions we repeat the simulations adopting now Ati = 1000 and At 2 = 10000. 
We observe in Fig. 6 that the return map has changed slightly when compared to Fig. 1. In particular, it no 
longer shows Ctot = 1 in the time series, but remains periodic, signaling that the system has shifted to a different 
cycle due to the perturbations. Apparently the periods observed in Figs 1 and 6 are the same. The distribution of 
periods and transient times are currently under investigation, results will be published elsewhere. Fig. 6 remains 
the same by increasing At 2 , which guarantees that the differences with respect to Fig. 1 are not a transient effect. 

According to these results, the role the noise plays, if allowed enough time to perturb the system significantly, 
is to drive the system from one attractor to a nearby one, which suggests that there is a family of periodic 
attractors which can be very close to one another in phase space. These effects, however, take place on a time 
scale which is longer [O(10 3 — 10 4 ) time steps] than the lifetime of the mice [C(10 2 — 10 3 ) time steps]. The small 
perturbations are therefore of little importance to the CA dynamics in the biologically relevant time scale, as 
suggested in previous work [9]. 

The behavior of the autocorrelation functions in Figs. 3 and 5 is reminiscent of what is observed in glassy 
systems. Here the CA approaches a periodic attractor, being thereafter deflected by the small perturbations. In 
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Figure 7: Cii{t + t w ;t w ) for different waiting times before small perturbations. From bottom to top, t w = 100, 5000, 
10000. 



glassy systems, temperature drives the system from one local minimum to another, preventing it from getting 
stuck in local minima of the potential energy surface. As the physical age of the system (characterized by the 
value of the waiting time t w ) grows, it finds itself exploring deeper and deeper regions of a rough potential energy 
surface, diffusing towards equilibrium [20]. Due to the roughness of the potential energy, equilibrium is only 
attained on very long time scales and relaxation is very slow. The older the system is the more it gets confined 
to a restricted region of phase space and the time scales for relaxation get longer and longer. Eventually, if we 
wait long enough, the system equilibrates and the dynamics becomes stationary, losing sensitivity to the waiting 
time. This picture of an aging physical system is reminiscent to the loss of plasticity for adaptation in a living 
organism as it gets older. In either case, one observes a strong dependence on the waiting time t w , evident when 
measuring two-time quantities like autocorrelation functions and responses. Results for the CA model are shown 
in Fig. 7, where we see the decay of autocorrelations for three systems with different ages or waiting times. Note 
that the horizontal axis is the time difference between the total time and the waiting time. The three curves 
should collapse in the case that the dynamical evolution is stationary. 



3.3 Large immunizations 

What is the role of the large perturbations on the dynamics of the system? From previous studies we would expect 
the large perturbations to accelerate the aging process: while the random small perturbations would change the 
route to the natural attractor of the system, the large ones would reduce the transient time, a conjecture that 
would explain the loss of plasticity of the older mice [9]. The protocol adopted is the one described in Section 2, 
with six-layer perturbations every 35 time steps, always at the same sites. 

In Fig. 8 we compare the results obtained for the unperturbed system and those of the system subjected 
only to large immunizations starting at t w (note that the densities are now also plotted as functions of t — t w ). 
Somewhat surprisingly, the decrease of the correlations for large perturbations is small, when compared to the 
case of the small perturbations (compare with Fig. 3). In hindsight, however, this can be understood because the 
large perturbations are always produced in the same region. For very long times the correlations will attain a 
stationary regime whose active populations contain at least part of the region involved in the immunizations [9]. In 
other words, the driving produced by the large periodic perturbations is of a completely different nature than that 
of the small perturbations. The large perturbations seem to play a selective role: the cycles that the system can 
reach are restricted to those that contain at least part of the populations incorporated during the immunization. 
According to this picture, the aging effects observed in Ref. [9] could be simply related to the exploration of the 
phase space: the older the system the less possibilities of choosing new cycles it would have. In Section 4 we 
discuss this issue further. 

In Fig. 9 we compare the time evolution of the densities B\ and B 2 and the autocorrelation functions for 
the system subjected only to large perturbations and for the system with both small perturbations and large 
immunizations. The results, as expected, confirm the dominance of the small perturbations, over the large ones, 
in driving the system faster to a different attractor. The increase of the densities around t ~ 10 4 for L = 50 
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corresponds to the same finite size effects occurring in Fig. 3, being associated to the active populations which 
have been incorporated into the network by means of the immunization protocol. 



4 Concluding Remarks 

We have analyzed the dynamics of a cellular automata model for the B-cell repertoire which has reproduced 
experimental results of immunizations on mice. Our analyses provide a broader context in which memory and 
plasticity take place, as discussed in Ref. [9]. 
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deviation). Inset: evolution for t > 100 shows that the HD attains a stationary regime for sufficiently long times 
(standard deviation not shown). 



Defining quantities analogous to the autocorrelation functions used in glassy systems, we have shown that 
the biologically relevant phenomena takes place in the transient regime of the model. The dependence of these 
functions on both t and t w (as opposed to the difference t — t w only) is reminiscent of the "aging" behavior observed 
in glassy systems, despite the fact that the underlying dynamics of both systems are controlled by completely 
different mechanisms. 

Starting from different initial conditions the deterministic dynamics takes the system towards a family of 
long-period attractors. When subjected to random small perturbations (Fig 4), the system is driven towards 
a new attractor of the family, revealing that most of the noise is assimilated. However, only part of the large 
perturbations is incorporated, due to the mechanisms of activation and suppression, leading to a saturation of the 
learning process [9]. 

From the biological point of view, the history of the mouse (sample) will be written by the different antigen 
presentations (random perturbations), starting from its "birth" (initial condition). Since the system is large but 
finite there is a maximum amount of information it may incorporate. The closer it is to its "destiny" attractor, 
the less information it is able to learn, since the deeper it already is in a given basin of attraction. Therefore the 
biological aging corroborated by the experimental results may be simply a consequence of this dynamical feature. 

From the results obtained up to now there are evidences that the purely deterministic dynamics is therefore 
non-ergodic. Despite their large number, however, the family of periodic attractors occupy only a fraction of 
the phase space. Evidences of this compression in phase space has been obtained in the following computer 
experiment: selecting randomly two initial configurations (with the same initial concentration), we measured the 
Hamming distance between them as a function of time as both systems evolve without any perturbation. In Fig. 10 
we show the evolution of the average HD for 100 pairs (error bars are standard deviations). During the first 10 
time steps the HD decreases very rapidly, but for t > 100 we observe a quasi-stationary regime, reflecting the slow 
driving to the attractors (see inset). Note that this behavior is somewhat similar to that of the autocorrelation 
functions for small t w . 

Still focusing on the evidence of phase space compression, in Fig. 11 we present the distribution of HD sampled 
from 500(500 — l)/2 pairs, for two consecutive time steps (t = 4000 and t = 4001). The distributions can be well 
described by a Gaussian, with a width that remains approximately constant (even for long times). Notice that the 
average value for t = 4001 is slightly larger than for t = 4000, reflecting the oscillatory behavior of the average HD 
as the pair of samples reaches their periodic attractors. It should be noted that, for large N, the Central Limit 
Theorem assures that randomly chosen initial configurations (with the same x) naturally give rise to a Gaussian 
distribution of the HD between any two of them, at time zero. Interestingly, the CA dynamics does not change 
the shape of the distribution, its only effect being to essentially shift the Gaussian towards lower mean values, on 
a long time scale. The spatio-temporal structure of the cycles, as well as their transients and basins of attraction, 
would be the object of further study and will be published elsewhere. 

This work was partially supported by CNPq, Faperj, Finep, Capes and Projeto Enxoval-UFPE. MC and RMZS 
acknowledge IF-UFF, where part of this work was done during their stay there. 
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